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SUPERSONIC FLOW CALCULATIONS FOR A CONE WITH AN ELLIPTIC FLARE"' 

by 

Harry Lehrhaupt 

A three-dimensional supersonic flow program, developed by the Grumman 
Corporation, has been used to calculate the flow field about a cone at zero 
angle of attack with an elliptical flare. A schematic diagram of the body 
is shown in Fig. (1). Calculations are carried out to two cone lengths, with 
the cone length taken to be unity. Free stream conditions of Mm = 8 
p = .00308308, and pm = 1 atmosphere were utilized in the calculations. 

Cross flow results in terms of the ratio W/v where v is the total velocity 
are shown in Figs. (2) to (5) for the meridian planes 0 = 20, 40, 60 and 80 
degrees, respectively. The pressure coefficient distributions 

C p = (P/Po, ' y/2 

on the body are shown in Figs. (6), (7) and (8) for the meridional planes 

9 = 30, 60 and 90, As a comparison in Figs, (6) to (8) C^ results for an 

2 

equivalent axisymmetric body utilizing Lomax's program have also been in^ 
eluded, An equivalent axisymmetric body is defined as that axisymmetric body 


* This work was supported by the National Aeronautics and Space 
Administration under NASA Grant NGR 33-016-131 , 

1) Three-rDimensional Near Characteristics program written by Dick Scheuing 
an employee of Grumman for his doctoral thesis, 

2) The program referred to here is the one described in the following report, 
Mamori Inouge, John V, Rakich and Howard Lomax, "A Description of Numerical 
Methods and Computer Programs for Axisymmetric Supersonic Flow Over Blunts 
Nosed and Flared Bodies," NASA TN-2970, August 1965, 
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which is tangent to the given body along a raeridinal plane. The distribu- 
tions for the equivalent bodies are seen to be in good agreement with the 
three-dimensional calculations. Such an agreement is to be expected at such 
high Mach numbers. 

In Figs. (9) and (10) the shock shape as computed by the 3-D calculations 
is compared with the axisymmetric calculations of the Lomax program for 
equivalent bodies at 0 = 60 and 9 = 90. The agreement between the shock shapes 

is even better than the agreement in the C distribution. 

P 

A user's manual for the three-dimensional program is included as an 
appendix to this memorandum. It should be noted here that the program, in 
its present version, is irrotational. As such the results obtained remain 
strictly valid in the region ahead of the first reflected characteristics 
from the points on the shock where the shock is no longer axisymmetric. For 
slowly varying geometries calculations could be extended outside of the above 
region, in particular under expansion. Special care must be taken, however, 
in adverse pressure situations where embedded shocks could be formed. 
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USER'S MANUAL 


Three Dimensional Near Characteristics 


The program has been broken down into the following subdivisions: 


A) GRUM 

B) Subroutine FIRST 

C) Subroutines FITl 

FIT2 

FIT3 


D) Subroutine CHAR 


E) FUNCTIONS VAL4, VAL 


F) Function HARDIF 

G) Subroutine NOMDIM 


H) Subroutines GRPH 
PLTS 


main program. 

reads inputs and initializes data. 

Spline fitting routines FITl is the most 
general. 

Calculates speed of sound squared and slopes 
of characteristics. 

written inorder to enable Scope to handle 
four dimensional arrays, 

3rd order, finite difference polynomial, 
nondimens ionalizes free stream velocities, 
pressures and densities, y = v /-/p»/p~ 

P " p /P°°; P = P/p*,- 

plotting routines W/v and CP are plotted 
as functions of X for pts. on the body. 
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INPUTS 


Inputs are broken down into three broad classifications RED, WHITE, BLUE. 


A) RED 

B) WHITE 

C) BLUE 


body configuration geometry, 
initial conditions. 

control parameters, step size tolerances, 
printout and plotting options. 


Coordinates used throughout are (X, r, 8). 



X 


Fig. A 
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RED INPUTS 


CARD 1 215, 17A4 


MIB 


MKB 


number of (r, 0) planes required to define 
the body. 

number of radii in defining body in each (r, 8) 
plane. 



CONFIG 


description of body being used (example: 
ELLIPTICAL CONE A/B 1. 394) 


CARD 2 (8F10. 5) 


XB(IB) IB= 1, MIB 
CARD 3 (8F10.0) 

TB (KB) KB= 1, MKB 


jjdlB/sl + 1 cards 

r i 

L J indicates greatest integer Example : 

[ 7 / 3." 1 = [3/5] = 0) 

stations where shape of body is defined, 

r 1 

L g j + 1 cards 

angles for radii in (r, 0) plane 0° <: 0 £ 90° 


or 0° <; 0 <; 180°. 
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CARD 4 
RB (IB, KB) 

CARD 5 
RBP(IB, KB) 


(8F10, 0) MIB*([*fT] + 1^) cards 

KB=1, MKB; IB=1, MIB - rsadii at each r-0 plane 

(8F10.0) 2 *(C^pJ + *) cards 

KB=1, MKBJ IB=1 and IB=MIB - slope ^ at first and last r-0 plane 
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WHITE INPUT 


CARD 6 
MJD 

HKD 

IVAXIS 

START 

CARD 7 

EMINF 

PINF 

RHQINF 

GAm 

ALPHA 

XI 


CAM I 
TB(KB) 

R§i(K8) 
RIXI (KB) 


(315, 16A4) 

number of points between body and shock 
(r) at initial profile, 
number of meridional points (g) data 
is specified at. 

velocities are cylindrical (1) or spherical 

( 2 ). 

alphanumeric label containing 64 characters, 


(6F12.0) (free stream flow conditions) 

free stream Mach number, 
free stream pressure, 
free stream density. 

C 

- y * ratio of specific heats , 

c v 

- angle of attack, 

- station at which we define initial velocity 
distribution, 

(6F1Q.Q) + l cards 

= | value of shock pt, (note at least two 

values are required; i.e,, MKB ^ -1). 

= r position of the shoek point. 

= slope of sheek |=, 
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CARD 9 


(4E15.0) MKD*MJD cards 


UI(JD,KD) 


VI(JD,KD) 
WI(JD,KD) 
PHI(JD, KD) 


x component of the velocity at pts. along 
initial profile, 
r component of the velocity. 

0 component of the velocity, 
helps locate where velocities are given 
between shock and the body 0 S PHI — ^ 
where PHI = 0 on the body and PHI = 1 
at shock 
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BLUE INPUTS 


CARD IQ (1015) 

MI 

MJ 

MK 

MAXM 

M\XN 

MODF 

IPT 

KPT 

I PUNCH 
NOPT 


number of steps we wish to march down- 
stream from the initial profile (Note: 
program terminates either when MI is 
exceeded or when x > XEND). 
number of pts. between body and shock 
(max is ?1 } 11 is coarse grid) (Note: 

MJD — MJ). 

number of meridional (0) planes (MK>MKD). 

=3 number of times we average slopes 
of characteristics. 

= 15 maximum number of iterations that, 
can be performed to find shock. 

(1) one used (affects iteration of shock pt. ) 

( 2 ) 

indicates r-0 plane at which we want total 
pressure calculated. 

(1) initial plane. 

(2) at each subsequent station, 
meridional plane at which total pressure 

drop at shock is computed 

(0) no punch or plot. 

(1) punch cards at last station and plots 

C and W/v versus X. 

P 

Print out occurs at every NOPT stations. 


9 



CARD 11 


(6F10.0) 


TOLEP 

.00001 tolerence for shock. 

CUNST 

0.8 (j/on Neuman stability constant.^ 

CTST 

0. 9 e stability factor 

CNX 

1.0 

CNNX 

0. 9 

XEND 

position where calculations are stopped 



THREE DIMENSIONAL NEAR CHARACTERISTIC PROGRAM 


Explanation of the Printed Output 


Complete Body Radius Matrix 
IB 
XB 
K 

RB 

RPB 

RPPB 

RPPPB 


Body cross-plane index number 
Axial location of a body cross-plane 
Meridional plane index number (0°-*180° or 
0° 90°) 


Body radial coordinate 

First derivative of r, 

b 

K - 22_) 

jx 

Second derivative (r^ 1 
Third derivative (r^ 11 


with respect 


2 

S r b 
3 


dx" 
_ r 


lx 5 


b_ ) 


to x 


Upstream Flow Conditions 
K 

Theta 

UU 

VU 

WU 


Meridional plane index number 
Meridional plane coordinate (0, in degrees) 
Axial component of the velocity upstream 
of the shock 

Radial component of the velocity upstream 
of the shock. 

Transverse component of the velocity 
upstream of the shock 
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New Data Surface Variables -Final Iteration 


I 

X2 

M 

K 

TH 

R2 

RX2 

RT2 

DR2 

DRT2 

RS2 

RSX2 

RST2 

J 

RP 

U2 

V2 

W? 

LAMBDAl 

LA.MBDA2 


New data surface index number . 

Axial coordinate of the new data surface* 
Number of iterations used in the calculation 
(M=l first order, M=2, 3 , 4 second order)* 
Meridional plane index number. 

Meridional plane coordinate (Q, in degrees)* 
Body radial coordinate (r^)* 

Body slope in the X direction ( ) * 

Body slope in the 0 direction ( •jg- ) • 

Radial distance between successive grid 

points between the body and the shock (6^)* 

Partial derivative of 6 with respect to 0 • 

Shock radial coordinate (r g > 

Shock slope in the X direction (|£) 

r dx s 

^ IT 

Shock slope in the 0 direction (^|-) ^ • 

Grid point index number in the radial 
direction (J=l on body). 

Local radial coordinate. 

Local axial velocity component. 

Local radial velocity component. 

Local transverse velocity component. 

Local slope of the first characteristic. 
Local slope of the second characteristic 
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New Data Surface Variables-Final Iteration (continued) 


MACH NO - Local Mach number 

CP - Local pressure coefficient * 


P/p 


CO 


y_ 

2 



1 


* Based on a constant total pressure drop across the shock for all points. 


13 





.01 X) 
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w/v (XIO -2 ) 












Fig. 5 
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Fig. 7 
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Fig. 10 
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3-D CALCULATION 
A-S CALCULATION 
BODY RADIUS 


X 


1.70 


1.90 


2.10 



APPENDIX A 


Listing of input data used to generate output for this report 
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780 

13 10 CIRCULAR CONE FLARE A/R-l 


0.0 

1.65 

1.0 

1.7 

1.1 

1.8 

1.2 

1.9 

1.3 

2.0 

1 .4 

1.5 

1.6 

0.0 

80.0 

10.0 

90.0 

20.0 

30.0 

40.0 

50.0 

60.0 

70.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.14482 

0.14577 

0.14485 

0.14580 

0.14493 

0.14506 

0.14522 

0.14539 

0.14555 

0.14568 

0.15798 

0.16178 

0.15810 

0.16191 

0.15843 

0.15894 

0.15957 

0.16025 

0.16090 

0.16143 

0.17115 

0.17969 

0.17140 

0.17998 

0.17211 

0.17323 

0.17463 

0.17617 

0.17764 

0.17887 

0.18431 

0.19947 

0.18473 

0.2 

0.18596 

0.18789 

0.19034 

0.19305 

0.19571 

0.19796 

0.19748 

0.22112 

0.19810 

0.2220 

0.19993 

0.20284 

0.20659 

0.21081 

0.21502 

0.21865 

0.21064 

0.24461 

0.21150 

0.24595 

0.21401 

0.21803 

0.22329 

0.22934 

0.23549 

0.24088 

0.21723 

0.25682 

0.21820 

0.25842 

0.22107 

0.22568 

0.23175 

0.23878 

0.24600 

0.25238 

0.22381 

0.26902 

0.22490 

0.27088 

0.22812 

0.23332 

0.24020 

0.24822 

0.25650 

0.26387 

0.23697 

0.29343 

0.23831 

0.29582 

0.24224 

0.24861 

0.25711 

0.26711 

0.27752 

0.28685 

0.25014 

0.31784 

0.25171 

0.32075 

0.25635 

0.26391 

0.27402 

0.28599 

0.29853 

0.30984 

0.26330 

0.34225 

0.26511 

0.34568 

0.27046 

0.27920 

0.29093 

0.30487 

0.31954 

0.33282 

0.13165 
0.1 3165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.13165 

0.24410 

0.13403 

0.24933 

0.14113 

0.15292 

0.16911 

0.18883 

0.21013 

0.22986 


21 2 1 7.5 DEGREE CONE FLARE-MACH 8-ALPHA 0 

8.0 2116.22 00308308 1.4 0.0 1.0 

0.0 0.189336 0.189336 90. C 0.189336 0.189336 

7662.5452484 1003.0345164 0.0000000 0.0000000 

7665.2918489 982.2874374 0.0000000 .0495957 

7667.9937256 962.3296744 0.0000000 .0992289 

7669.6591866 943.0656566 0.0 0.1489004 

7673.2959001 924.4107210 0.0000000 .1986111 

7675.9110928 906.2887550 0.0000000 .2483619 

7678.5117211 888.6302815 0.0000000 .2981535 

7681.1046278 871.3708484 0.0000000 .3479867 

7683.6966936 854.4496187 0.0000000 .3978626 

7686.2946918 837.8080570 0.0000000 .4477818 

7688.9069613 821.3886714 0.0000000 .4977452 

7691.5406049 805.1336870 0.0000000 .5477537 

7694.2047323 788.9835793 0.0000000 .5978081 

7696.9092701 772.8753366 0.0000000 .6479093 

7699.6656742 756.7402754 0.0000000 .6980581 

7702.4874993 740.5011376 0.0000000 .7482554 

7705.3912149 724.0680319 0.0000000 .7985021 

T*«e. 3974240 707.3324669 0.0000000 .8487991 

7711.5327652 690.1581140 0.0000000 .8991471 
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7714.8330479 
7718.3487679 
7662.5452484 
7665.2918489 
7667.9937256 
7669.6591866 
7673.2959001 
7675.9110928 
7678.5117211 
7681.1046278 
7683.6966936 
7686.2949918 
7688.9069613 
7691 . 5406049 
7694.2047223 
7696.9092701 
7699 . 6656^42 
7702.4874093 
7705 . 3912149 
7708.3974240 
7711.5327652 
7714.833079 
7718.3487679 
150 21 10 

0.00001 0 . 


672.3656533 0 

653.7061513 0 

1003.0345164 0 

982.2874374 0 

962.3296744 0 

943.0656566 0.0 

924.4107210 0 

906.2887550 0 

888.6302815 0 

871.3708484 0 

854.4496187 0 

837.8080570 0 

821.3886714 0 

805.1336870 0 

788.9835793 0 

772.8753366 0 

756 . 7402 7 54 0 

740.5011376 0 

724.0680319 0 

707.3324669 0 

690.1581140 0 

672.3656533 0 

653.7061513 0 

3 15 1 2 

8 0.9 1.0 


0000000 

.9495471 

0000000 

1.0000000 

0000000 

0.0000000 

0000000 

.0495957 

0000000 

.0992289 

0 . 

1489004 

ooooooo 

.19861-11 

0000000 

.2483619 

ooooooo 

.2981535 

ooooooo 

.3479867 

ooooooo 

.3978626 

ooooooo 

.4477818 

ooooooo 

.4977452 

ooooooo 

.5477537 

ooooooo 

.5978081 

ooooooo 

.6479093 

ooooooo 

.6980581 

ooooooo 

.7482554 

ooooooo 

.7985021 

ooooooo 

.8487991 

ooooooo 

.8991471 

ooooooo 

o 949 54 7 1 

ooooooo 

1.0000000 

1 1 

10 

0.9 

1.36 
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APPENDIX B 


Program listing 
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/ 


Y8 2 5007. PI .T1300.CM13700n.LEHRHAUPT(GRUM-PLOT ) 

PUN ( S ) 

LGO ( L077000) 

780 

8L0CK DATA 

COMMON/ T I TLE/TITL1(3).TITL2(3) .LABLEU ) /PLTS/XORI .YOR I .SX.SY* J 
1L* IT PN .VALUE. H 

DATA TI TL1 .T ITL2 ,XORI .YORI ,SX .SY. J.L.ITPN/30H X 

1 * 30H ,0,5.0. 5. 15.0.9 

2,3,8/ » L ABLE , VALUE /3 OH CROSS FLOW THETA a .0.12.56/ 

3 0.2/ 

END 
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rv r\ n n 


PRECEDING PAGESBLANK NOT FILMED 


DS = DT 

CALL F I T 3 ( MK * 1 ) 
DO 2314 K=1,MK 
2314 RT2 (K)=FP(K) 


FIRST ORDER SOLUTION (M=l) 
M= 1 

DO 3214 K = 1 , MK 
DRST=CVNST*DR1 < K ) 


SHOCK POINT (M=l) (3000) 


N= 1 
I N 1 = 1 
I N 2 = 1 

RSTR=RST1(K)/RS1(K) 

RSX ( 2 ) = R S X 1 ( K ) 

RC=RS1 ( K ) -DR1 ( K ) 

EKAY1=RC+DX*EL1(MJM1»K,1) „ .... 

EKAY2*- 1 , / ( DR1 ( K ) +DX* ( EL 1 ( MJ .K ♦ 1 l-ELl (M JMl »K * 1 ) > 1 

2001 Rd=RS1(K)+.5*(RSX1{K)+RSX(2)) *DX 
EK A Y = ( RP-EKAY 1 ) *EKAY2 
RA ( 1 )=RC+EKAY*DR1(K) 


2 002 FM1?NN)=VAl!fi. 1*MJM1.K,1»NN)+EICAY*(VAL(F1.1*MJ.IC.1*NN) VAUFl.l. 

^ CALL CHAR <FA(1.1)*FA(1*2)»FA(1.3)*0**A*ELA(1»1)»ELA(1»2)) 

GA( 1 ) = ( 2.*FAC 1.3)*(FA<1 *1)*FA( 1 ,4)+FA( 1»2)*FA(1 *5) ) (A-FA< 1,3) 2) 

j ( 1,6) -FA (1»2)*<A+FA(1 .31**2 > ) / ( RA (1) *< A-F A I 1 , 1 > ** 2 ) ) 

;.,GA ( 1 ) = FA ( 1 » 1 ) +ELA( 1 » 2 )*FA( 1 » 2 >+GA ( 1 )*DX 

to I A = ( F A ( 1 » 4 ) +E LA ( 1 ♦ 1 ) * ( FA ( 1 » 5 1 -FA (1.3) > ) /RA < 1 ) 

tfW A=FA( 1.3) +W1A*DX 

W2CMJ»K)*WWA 

3 003 ENRS=1./SQRT(1.+RSX<2)**2+RSTR**2) 

GO TO (3004.3005.3004) .INI 

30 04 VNU=ENRS*(RSX(2)*UU<K)-VU(K)+RSTR*WU<K) 1 


VNU2 = VNU**2 
t-PS = G2+G3/VNU2 

3005 V2 ( MJ »K ) =VU ( K > + < 1 . -EPS ) *ENRS*VNU 
J 2 ( M. J . K ) =GGA (1)-ELA(1»2)*V2(MJ.K) 

VD2 = U2 ( MJ »K ) **2+V2 (MJ ,K ) **2+W2 ( MJ »K ) **2 

VND=SQRT (VD2+VNU2-VINF2 ) 


FPC=VND/VNU 

OEPEP(2)=(EPC-EPS)/EPS 

IF (TOLEP-ABS(DEPEP( 2 ) ) ) 3006.3011.3011 

3006 GO TO (3007,3008) .TN2 

3007 I N 1 = 2 
IN2 = 2 

EPS=.5*(EPS+EPC) 

VNU2 = G3 / ( EPS-G2 ) 

VNU=SQRT (VNU2 ) 

RSX(2)= l (UU(K ) )*(VU(K)-WU(K)*RSTR) + VNU*SORT( (VINF2~VNU2)*(1. + RSTR**2 
1 ) — ( VU ( K ) *R5TR+WU ( K ) )**2 ) ) / <UU(K)**2 - VNU2> 


GO TO 3009 


3008 


I N 1 = 3 


TFMPOR = RSX< 2 ) 
RSX(2)=RSX(2>-(RSX(1 


) -RSX ( 2 ) )*DEPEP<2) /(DEPEP(1 )-DEPEP(2) 1 


RSX ( 1 )=TEMPOR 
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3009 DEPEPf 1 )=DEPEP(2 ) 

N = N+ 1 

IF (N-MAXN) 3010*9001 * 900 1 

3010 GO TO ( 3001 ,3003 ), MODE 

3 011 RS ?(<)= RSI ( K ) + .5*(RSXl (K )+RSX ( 2 )) *DX 
95X2 ( K ) =RSX ( 2 ) 

CALL CHAR (U2 (MJ,K ) ,V2(MJ,K),W2(MJ,K)*0.,A*El-2(MJ,K,l),EL2(MJ,K,2)) 
EL IB < MJ,K ) = . 5*( EL A { 1 , 1 ) +EL2 ( MJ , K » 1 ) ) 

E L 2 B ( M J , < ) = E L 2 ( M J » K , 2 ) 

DR2(K)=(RS2(K)-R2(K) )/MJMl 

c 

C, BODY POINT ( M= 1 ) (3100) 

RP = R2 ( K ) 

RC=R1 ( K ) 

EKAY= (RP-RC-DX*EL1 ( 1 , K , 2 ) ) /(DR1 (K )+DX* (ELI ( 2»K»2 ) -EL 1 ( 1 * K * 2 ) > ) 
RA(2)=RC+EKAY*DR1(K) 

DO 3101 NN= 1,6 

3101 FA ( 2 » NN ) = VAL ( F 1 » 1 » 1 » K » 1 »NN ) +EKAY* ( VAL ( F 1 » 1 » 2 » K » 1 »NN ) - VAL ( F 1 » 1 » 1 » K * 
1 1 , N N ) ) 

CALL CHAR ( F A ( 2 » 1 ) , FA ( 2 * 2 ) *FA( 2 * 3 )» 0 .,A*ELA( 2 * 1 ),ELA( 2 » 2 ) ) 
GA{ 2 )=( 2 o*FA( 2 , 3 )*(FA( 2 ,l}*FA< 2 » 4 )+FA( 2 » 2 )*FA( 2 , 5 ))-(A-FA( 2 » 3 )** 2 ) 
1 *F A ( 2 » 6 ) ~FA ( 2 » 2 )*(A+FA( 2 » 3 )** 2 ) ) / ( RA ( 2 ) * ( A-FA ( 2 , 1 ) ** 2 ) ) 

GGA ( 2 ) -FA ( 2 , 1 ) +ELA ( 2,1) *FA (2,2) +GA ( 2 >*DX 
W2B=(FA(2»4)+ELA(2,2)*(FA(2,5 )-FA (2*3) ) ) /RA ( 2 ) 

WWB=FA(2»3)+W2B*DX 

V2 ( 1 * K ) = ( GGA ( 2 )*RX2 ( < )4-WWB*RT2 (K.) /R2 ( K ) ) / ( 1 • + RX2 ( K ) *ELA (2,1)) 

U2 ( 1 , K ) =GGA ( 2 ) -ELA ( 2 , 1 )*V2 ( 1 ,K) 

W2 ( 1 , K ) = WW8 

CALL CHAR(U2 ( 1 , K ) , V2 ( 1 , K ) , W2 ( 1 » < ) , 0 . , A » EL2 ( 1 , K , 1 ) , EL 2 ( 1 , < , 2 ) ) 

EL 1 B ( 1 ,K)=EL2(1»K,1) 

EL23(i,K) = «5*(ELA(2»2)+EL2(l»K»2) ) 

r 

C FIELD POINTS (M=l ) (3200) 

DO 3213 J = ? »MJMl 
RP = RP+DR2 ( K ) 

L = 0 

JJ=J-1 

3201 L = L + 1 

3202 RC=R1 (K )+( JJ-1 )*DR1 (X ) 

IF ( RP-RC-DX*EL1 ( JJ,K ,L ) )3203,3204,3205 

3203 JJ=JJ-1 

GO TO 3202 

3204 EKA Y = 0 • 0 
GO TO 3208 

3205 IF ( RC+DR1 ( K )+DX*ELl ( JJ+1 ,K ,L )-RP ) 3206,3206,3207 

3206 JJ=JJ+1 

GO TO 3202 

32 0 7 EKAY=(RP-RC-DX*EL1 ( JJ,K,L) ) / ( DR 1 ( K ) +DX* ( EL 1 ( J J+l ♦ K , L ) -EL 1 ( J J , K , L ) ) 

1 ) 

3208 RA(L)=RC+EKAY*DR1(K) 

DO 3209 NN= 1,6 

32 09 FA(L ,NN)=VAL(F1»1 , J J , K , 1 , NN ) +E< A Y* ( VAL { FI , 1 , J J+l , K , 1 , NN ) -VAL ( F 1 , 1 , 
1JJ»K , 1 ,NN ) ) 

CALL CHAR(FA(L,1) , FA ( L , 2 ) , FA ( L , 3 ) , 0 . , A , ELA ( L , 1 ) , ELA ( L , 2 ) ) 
GA(L)=(2.*FA(L,3)*(FA(L,1)*FA<L,4)+FA(L,2)*FA(L,5) ) - { A-FA ( L , 3 ) **2 ) 
1*FA'(L*6)-FA(L.2)*<A+FA(L.3)**2> ) / ( RA ( L ) * ( A-FA ( L , 1 ) **2 ) ) 

GGA (L ) =FA (L , 1 ) +ELA ( L , 3-L ) *FA ( L , 2 ) +GA(L)*DX 
GO TO ( 3201 ,3210) »L 
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3210 IF (DRST~RA(2)+RA(1 ) ) 3 2 1 1 ♦ 32 12 . 3 2 1 2 

3211 CNX=CNX*CNNX 
GO TO 2305 

3212 W1A- (f A ( 1*4 )+ELA( 1 *1 )*(FA ( 1 *5 5 -FA (1.3) ) ) /R A ( 1 ) 

WW A = F A ( 1*3) +W1A*DX 

V2 ( J.K) = (GGA(2)-GGA( 1 ) ) / ( EL A ( 2 * 1 ) -EL A ( 1.2) ) 

U2 < J .K) =GGA ( 1 5-ELAC 1 .2 )*V2 ( J.K) 

W2(J.K)=WWA 

CALL CHAR (U2 (J.K) ♦ V2 < J ♦ it ) . W2 ( J . K ) ♦ 0 . ♦ A * EL2 ( J . K . 1 ) . EL 2 < J ♦ K ♦ 2 ) ) 
EL1BI J.K )s.5*( ELA ( 1 .1 ) +EL2 ( J *K ♦ 1 ) ) 

3213 EL28( J ♦ K ) = . 5* ( ELA ( 2 ♦ 2 > +EL2 ( J.K. 2) ) 

3214 CONTINUE 

DO 5209 M=2.MAXM 

C COMPLETE NEW data SURFACE MATRIX (4000) 

DS = DT 

DO 4001 K = 1 »MK 

4001 F(K)=DR2(K) 

CALL F I T3 ( MK » 1 ) 

DO 4002 K= 1 .MK 
DRT2 ( K ) =FP ( K ) 

4002 RST2(K)=RT2(K)+MJM1*DRT2(K) 

DO 4104 K=1.MK 

DS = DR2 ( K ) 

DO 4103 NN = 1 ♦ 5 
DO 4101 J=1,MJ 

4101 F< J)=VAL(F2.2.J.K.1.NN) 

CALL FIT1(MJ.3,3,0..0. ) 

DO 4102 J=1.MJ 

4102 CALL VAL4(F2.2,J*K,3.NN»FP< J) ) 

4103 CONTINUE 

4104 CONTINUE 

DO 4206 J = 1 » M J 
DS = DT 

DO 4205 NN = 1 .5 
KEND= 1 

IF (3-NN) 4202,4201.4202 

4201 KEND=2 

4202 DO 4203 K=1,MK 

4203 F(K)=VAL(F2.2.J.K,1.NN) 

CALL F I T3 ( MK » KEND ) 

DO 4204 K=1 ,MK 

4204 CALL V AL4 ( F2 ♦ 2 , J . K , 2 » NN , FP ( K ) -VAL ( F 2 » 2 .J.K.3»NN)*(RT2(K)+(J-1)* 
1DRT 2 ( K ) ) ) 

4205 CONTINUE 

4206 CONTINUE 
MO=M-l 

C 

C SECOND ORDER SOLUTION (M=2»3....) 

DO 5208 K = 1 . MK 
DRST=CVNST*DR1(K) 

c 

c SHOCK POINT ( M = 2 * 3 . • • • ) (5000) 

N= 1 
INI = 1 
I N2 = 1 

RPO=RS2 ( K ) 

RSTR=RST2(K)/RS2(K) 
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KOX ( ^ ) = KSX2 ( K I 

A = A0-G1* (U2 ( MJ ,K ) **2+V2 <MJ,K)##2+W2(MJ,K)**2) 

GP= ( 2.*W2 (MJ,K )* (U2 (MJ»K ) *UT2 ( M J ♦ K ) +V2 ( MJ , K > *VT 2 ( MJ » K ) ) 
l-( A-W2 (MJ»K ) **2 ) *WT2 < MJ , K) - V2 ( MJ , K ) * ( A + W2 ( MJ ,K ) **2 ) ) / (RPO* 

2( A-U2 (MJ,K) **2 ) ) 

W1P=(UT2(MJ»K)+EL2(MJ »<»!)*{ VT2(MJ.K)-W2(MJ»K) ) )/RPO 
5001 RP = RS1(K)+.5*(RSX1(!<)+RSX(2) )*DX 
RA ( 1 ) =RP— V AL ( F2 * 2 ,M J ♦ K * 1 , A ) *DX 
CAPA ( 1 ) = -VA L ( F 2 , 2 , ^ J * K * 2 ♦ A ) *DX/RA ( 1 ) 

DR = RA ( 1 ) -RSI (K)+DRl ( K ) 

DO 5002 NN= 1*6 

5 0 02 C A( 1 ,NN ) = VAL <F1 , 1 *MJMl * K , 1 , NN ) +DR * ( VAL ( F 1 , 1 * M JM 1 * < , 2 t NN ) + , 5 *DR* 
1(VAL(F1»1 *M JM1 » K » 3 » NN ) +DR* VAL(F1*1 » MJM 1 ♦ K , A , NN > /3 . 0 ) ) 

CALL CHAR (FA(1,1),FA(1»2)»FA(1»3) »CAPA (1),A»ELA(1,1) t FLA < 1 , 2 > ) 

DO 5003 NN = 4 *6 

FAR =VaL (FI* 1 , MJM 1 , < , 2 , NN-3 ) +DR* ( VAL (Fit 1 ,MJMl »K,3»NN-3)+ 

10. 5*DR*VAL( FI , 1 t MJM ] , K , 4 * NN-3 ) ) 

5 0 03 FA ( 1 > NN ) =FA ( 1 , NN ) - VAL (F2 t 2 t MJ ,K t 2 .4 5 *DX*FAP. 

GA ( 1 ) = ( 2 • *FA ( 1 ,3 ) *FA<1 1 1 )*FA ( ! ) + ( 2.*FA< It 2 ) *FA ( ] ,3 )+( A-FA( 1 ,? ) 

1**2 ) *CAPA( 1 ) )*FA ( 1 ,5 ) -( A-FA ( 1 ,3 )**2 )*FA ( 1 ,6 ) ~ ( FA ( 1 * 2 ) *FA ( 1 , 3 ) + 
2(A-FA(1i3)**2)*CAPa( 1) ) *FA ( 1 t 3 ) -A*FA ( 1 t ?) ) / ( R A ( 1 ; * ( A~FA ( 1 , 1 ) *«p ) ) 
GBAR ( 1 ) = .5* (GA ( 1 ) +GP) 

ELBAR( 1 ) =.5* (FLA ( 1 ,2 ) +EL2 (MJ.K t? ) ) 

GG A ( 1 ) = F A ( 1 1 1 ) •‘•ELBAR ( 1 ) *F A ( 1 t 2 5 +GRAR ( 1 ) *DX 
W1A= ( FA ( 1 ,4 )+ELA ( 1 ,1 ) * (FA ( 1 »5 J-^A ( 1 ,3 >) ) /RA( 1 ) 

WlBAR=. 5* (W1A+W1P ) 

WW A = F A ( 1 t 3 ) +.5*CAPA ( 1 ) *FA { ! ,2 )+Wl RAR*DX 
5 0 04 ENRS = 1 . /SORT ( 1 ,+PSX ( 2 ) **2+RSTR**2 ) 

C-0 TO (5005.5006,5005 ) t INI 
5005 VNU = ENRS*(RSX( 2)*UU(K ) -VU ( K )+RSTR*WU(K ) ) 

VNU2 = VNU**2 
EPS=G2+G3/VNU2 

5 006 V2 ( M J ♦ K ) =VU(K ) + ( 1 .-EPS )*ENRS*VNU 
U2 ( M J , K ) =GOA ( 1 ) -ELRAR ( 1 )*V2 (MJ,K ) 

W2(MJ,K) =WWA— ,5*CAPA( 1 )*V2(MJ,K) 

VD2=U2 ( MJ,K )**2+V2 (MJ ,K ! **2+W2 ( MJ ,K ) **2 

VND = SORT ( VD2+VNU2-VINF2 J 

EPC=VND/VNU 

DEPEP(2 )=(EPC-EPS)/FPS 

IF ( T OLEP-ABS ( DEPEP ( 2 ) M 5007,5 012,5012 

5007 GO TO ( 5008 ,5009) , IN2 

5008 I N 1 =2 
IN2 = 2 

EPS=.5*(EPS+EPC) 

VNU2 = G3/(EPS-G2 ) 

VNU=SORT ( VNU2 ) 

RSX ( 1 ) = RSX ( 2 ) 

RSX ( 2 ) = (UU(K )*{VU(K ) — WU ( K )*RSTR )+VNU*SQRT ( ( V I NF2~VNU2 ) * ( 1 .+RSTR**2 
1 ) - ( VU ( K ) *RSTR+WU ( K ) ) **2 ) ) / (UU( K ) **2~VNU2 ) 

GO TO 5010 

5009 I N 1 = 3 
TEMPOR=RSX( 2 ) 

RSX ( 2 ) =RSX ( 2 )-(RSX ( 1 ) -RSX ( 2 ) )*DEPEP ( 2 ) /(DEPEP ( 1 )~DEP EP( 2 ) ) 

RSX ( 1 ) = TEMPOR 

5010 DEPEP(1)=DEPEP(2) 

N = N + 1 

IF (N-MAXN) 5011,9001,9001 

5011 GO TO ( 5001 , 5004 ), MODE 
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5012 RS2 (K ) = RS1 ( K ) + . 5* (RSX1 (K )+RSX { 2 ) ) *DX 
RSX2 ( K ) =RSX ( 2 ) 

CALL CHAR t U2 ( MJ » K ) , V2 (M J , K ) , W2 ( MJ , K ) *0 * , A , FX2 ( M J »K * 1 ) ,EL2 ( M J , K , 2 ) ) 
EL IB (MJ*K ) = .5* (EL A ( 1 ♦ 1 )+FL2 (MJ»K» 1 ) ) 

EL2B ( MJ * K )=EL2 ( MJ ,K » 2 ) 

DR0=DR2(K) 

DR2(K)=(RS2(K)-R2(K) )/MJMl 
C 

C BODY POINT ( M = 2 »3 * • • • ) (5100) 

RP0=R2(K) 

A=A0-G1* (U2 ( 1 *K )**2+V2 ( 1 ,K )#*24-W 2 ( 1 »K) **2) 

GP= (2.*W2 ( 1 *K) * (U2 (1 *K )*UT2 ( 1 ,K )+V2 ( 1 ♦ K ) *VT 2 ( 1 , K )) - ( A-W2 ( 1 * K ) **2 ) * 
1WT2(1«K}-V2(1*K)*(A+W2(1*K}**2) ) / ( RPO* ( A-U2 * 1 ,K)**2) ) 

W2P=(UT2 ( 1 *K )+EL2 {1*K,25*(VT2(1*K)-W2(1,K)!) /RPO 
RP=R2 ( K ) 

RA ( 2 ) =RP- VAL ( F 2 » 2 » 1 *K , 1 ,5 ) *DX 
CAPA ( 2 ) =-VAL ( F2 * 2 » 1,K»2»5)*DX/RA(2) 

DR=RA(2 ) -R 1 <K) 

DO 5101 NN=1,6 

5101 FA(2*NN)=HARDIF(VAL(F1 *1*1 ,K,1,NN) *VAL { FI * 1 » 1 ,K *2 ,NN ) , 

1 VAL ( FI » 1 » 1 » K » 3 »NN ) , VAL (F 1,1 * 1 , K , 4 * NN ) * DR ) 

CALL CHAR ( FA < 2 , 1 ) , FA ( 2 * 2 ) , FA ( 2 , 3 ) , CAP A ( 2 ) , A , EL A ( 2 * 1 ) , EL A ( 2 , 2 ) ) 

DO 5102 NN=4,6 

FAR=VAL (FI , 1 , 1 ,K ,2 ,NN-3 )-DR*( VAL ( FI » 1, 1 ,K,3 »NN-3 )+0. 5*DR* 
lVAL(Fl,l,l,K,4,NN-3) ) 

5102 FA ( 2 »NN ) =FA ( 2 , NN ) -VAL (F2,2,1,K,2»5 )*DX*FAR 
GA(2)=(2.*FA(2»3)*FA{2,l)*FA(2»4)+(2.*FA(2*2)*FA(2»3)+(A-rA(2»3) 

1**2 )*LAPA(2) )*FA(2,5 ) - ( A-FA ( 2 » 3 ) **2 ) *FA ( 2 , 6 ) - ( FA ( 2 » 2 ) *FA ( 2 * 3 ) + 

2 ( A-FA ( 2 *3 )**2 ) *CAPA ( 2 ) )*FA( 2 »3 ) -A*F A ( 2 » 2 > ) / ( R A ( 2 ) * ( A-FA ( 2 » 1 ) **2 ) ) 
G8AR(2)=.5*(GA(2)+GP) 

ELBAR(2)=. C *(ELA(2,1 )+EL2( 1»K,1 ) ) 

GGA ( 2 )=FA( 2 ,1 )+ELBaR ( 2 ) *FA ( 2 » 2 ) +GBAR ( 2 )*DX 
W2B=(FA(2»4)+ELA(2,2)*(FA(2,5)-FA(2,3) ) )/RA(?) 

W2BAR=«5* ( W28+W2P ) 

WWB=FA( 2 ,3 )+.5*CAPA( 2 ) *FA ( 2 ,2 ) +W28AR*DX 

V2 ( 1 »K ) = (GGA ( 2 ) *RX2 ( K ) +WWB*RT 2 ( K ) /R 2 ( K ) > / ( 1 »+ELBAR ( 2 ) *RX 2 ( < ) + . 5* 
1CAPA(2)*RT2 (K) /R2 ( K ) ) 

U2 ( 1 » K ) =GGA (2)-ELBAR(2)*V2(l»K) 

W 2 C 1 * K ) =WWB-« 5*CAPA(2)*V2(1 ,K) 

CALL CHAR ( U2 ( 1 , K ) , V2 ( 1 » K ) , W2 ( 1 » K ) , 0 . , A » EL2 ( 1 , K , 1 ) , EL2 ( 1 , < , 2 ) ) 

EL1BI 1,K)*EL2< 1,K,1 ) 

EL2B(1,K)=.5*(£LA(?,2)+EL2(1,K»2) ) 
r 

C FIELD POINTS (M=2»3*...) (5200) 

DO 5207 J=2,MJM1 

IF (DRST-(EL2B( J,K)-EL1B( J,<) )*DX ) 5201,5202,520? 

5201 CNX=CNX*CNNX 
GO TO 2305 

5202 RPO=RPO+DRO 

A=A0-G1*(U2 ( J»K)**?+V2( J,K)**?+W?( J»K) **2) 

GP=(2.*W2( J»K)*(U2(J,K)*UT2( J,K)+V2( J»K)*VT2( J,K) ) - ( A-W2 ( J » K ) **2 ) * 
1WT2( J,K)-V2( J,K)*(A+W2( J,K)**2) )/ (RPO* (A-U2 ( J *K ) **2 ) ) 

W1P=(UT2( J,K)+EL2( J,K,1 )*(VT2(J»K)-W2(J»K>) ) /RPO 
RP=RP+DR2 ( K ) 

L = 0 

5203 L=L+ 1 

RA(L)=RP-VAL(F2,2,J,K,1 ,L+3)*DX 
CAPA(L)=-VAL(F2»2,J,K,2,L+3)*DX/RA(L) 
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n n 


JJ=1+INT( (RA(L)-Rl(K) )/DRl(K) ) 

DR = RA(.L )-Rl (K)-( JJ-1 ) *DRl ( K ) 

DO 5204 NN= 1 .6 

5204 FA(L.NN)=H/RDIF(VAL(F1 . 1 .JJ.K.l .NN) » VA L ( F 1 * 1 . J J ,K , 2 » NN ) ♦ 
1VAL(F1,1,JJ.K,3»NN) , VAL ( F 1 , 1 , J J ,K , 4 . NN ) * DR ) 

CALL CHAR (FA (L*l ).FA(L,2).FA(L*3) .CAPA (L) . A » ELA ( L , 1) . EL A ( L ♦ 2 I ) 

DO 5205 NN=4,6 

FAR=HARDIF( VAL (FI , 1 ,K ,2 ,NN~3 ) ,VAL (FI ♦ 1 * J J , < , 3 * NN- 3 ) , 

1VALIF1.1. JJ .K.4.NN-3) .0.0. DR) 

5205 FA(L,NN)=FA(L,NN)-VAL(F2,2 .J.K.2.L + 3 >*DX*FAR 

GA (L)=(2.*FA(L,3)*FA(L.1)*FA(L,4)-M2»*FA(L.2)*FA{L.3 )+( A-FA(L.3) 
1**2 )*CAPA(L) )*FA(L,5)-(A-FA(L.3)**2)*FA(L.6)“(FA(L»2 )*FA(L»3)+ 
2(A-FA(L,3)**2)*CAPA(L) ) *FA ( L . 3 ) - A*F A ( L , 2 ) ) / ( R A ( L ) * ( A-FA ( L » 1 ) **2 ) ) 
GBAR(L)=.5*(GA(L)+GP) 

FLBAR(L)=.5*(ELA(L,3-L ) +EL 2 ( J ♦ K ♦ 3-L ) ) 

GGA(L)=FA(L»1 )+£LBAR(L)*FA(L.2)+GBAR(L)*DX 
GO TO ( 5203.5206 ) ,L 

5 2 06 W1A= (FA( 1 ,4)+ELA( 1 ,1 ) *(FA ( 1 .5 )-FA (1.3) ) ) / RA ( 1 ) 

W 1 BAR= .5* ( W1 A+WlP ) 

WWA = FA ( 1.3) + .5*CAPA(1 ) *FA ( 1 .2 )+WlBAR*DX 
V2 ( J .K ) = ( GGA ( 2 ) -GGA ( 1 ) ) / ( EL BAR ( 2 ) -ELBAR ( 1 ) ) 

U2< J.X)=GGA( 1 ) -ELBAR ( 1 )*V2( J.K) 

W2 ( J *K ) = WWA- • 5 *CAP A ( 1 )*V2( J.K) 

CALL CHAR(U2 (J»K).V2(J*K).W2(J*K).0..A*EL2(J.K.1).EL2(J.K*2)) 
EL1B(J,K)=.5*(ELA(1,1 )+EL2( J.K.l) ) 

5207 EL2B( J,K)=.5*( ELA(2»2 )+EL2( J.K. 2) ) 

5208 CONTINUE 

5209 CONTINUE 

DELTA THETA STABILITY AND PRESSURE COEFFICIENT (6000) 

MO=M- 1 

GO TO (6001 .6002*6004 ) *IPT 

6001 IPT =3 

VNU2=(RSX1(KPT )*UU(KPT ) -VU ( <PT ) +RST 1 ( KPT ) *WU ( KPT ) /RSI (KPT ) )**2/ 

1( 1.+RSX1 (KPT )**2+(RSTl (KPT } /RSI (KPT) )**2) 

GO TO 6003 

6002 VNU2= (RSX2(KPT)*UU(KPT ) -VU( KPT ) +RST2 ( KPT ) *WU ( KPT ) /RS2 (KPT) )**2/ 

1( l.+RSX2(KPT )**2+(RST2 (KPT) /R S2 ( KPT ) ) **2 ) 

6003 PTRAT=(G5*VNU2-G2 )/( ( G2+G3/VNU2 ) * ( G5*VNU2-G2 ) ) **G4 

6004 NPRT=1 

I F ( I JL.NE.I JL/NOPT*NOPT )NPRT=»2 
I JL = I JL + 1 

GO TO ( 604 .607 ) , NPRT 
604 WRITE (IW. 20019) I.X2.M0 

20019 FORMAT ( 1H1 . 30X » 53HNEW DATA SURFACE VARIABLES - FINAL ITERATIO 
IN (1 = 15. 6H ♦ X2 = F9.5 *5H . M=I5»1H)/) 

607 DO 6007 K=1.MK 

GO TO (608*609) .NPRT 

608 WRITE (IW, 20020) K , T ( K ) , R 2 ( K ) , RX2 ( K ) . R T2 ( K ) , DR2 ( K ) ,DR T2 ( K ) , RS2 ( K ) , 
1RSX2 (K) , RST2 (K ) 

20020 FORMAT {///3H K=I4,4H TH*F9.4,4H R2=F10.6,5H RX2=F8.6,5H RT2=F9.6. 
15H DR 2 = F8 . 6 » 6H DRT2 = F9.6,5H RS2 = Fl0.6»6H RSX2=F8.6,6H RST2 = F9.6// 

2 14X , 1HJ.6X.2HRP.12X.2HU2 , 12X ,2HV2 , 1 2X . 2HW2 , 1 OX , 7HLAMBDA 1 .7X, 
37HLAMBDA2 ,7X .7HMACH NO.9X.2HCP/) 

609 RP=R2 ( K ) -DR2 ( K ) 

DO 6006 J = 1 » M J 
RP=RP+DR2 (K ) 

VP2=U2( J.K)**2+V2( J,K)**2+W2< J,K)**2 
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A=A0-G1*VP2 

EM2 =VP2 / A 

EM( J,K) «=SQRT (EM2 ) 

BETA=SQRT(EM2-1. ) 

TANALF=ABS( W2 ( J*K) ) /SORT (U2LL&& )##2+V2 ( J,K)**2 ) 

ELAMB=( BETA »T ANAL F+l . ) / ( 8ETA-TANALF ) 

TEST=DX/OT-CTST* (RP*U2 ( J 9 x )-OX^V2( J»K ) ) / ( EL AMB*SQRT ( U2 ( J * K ) **2+ 
1V2( J,K)**2) ) 

IF ( T EST ) 6005*6005*9002 

6005 CP ( J »K ) = { PTRAT*A*#G4-CCP1 ) *CCP2 
GO T0(6006*6007) *NPRT 

6006 WRITE (IW, 20021 ) J ,RP *U2 ( J * K ) , V2 ( J * K ) • W2 ( J . K ) , EL2 f J » K * 1 ) , 

1EL2( J»K»2) »EM( J.K) ,CP(J»K) 

20021 FORMAT ( 1 OX , 1 5 , 8 F14 • 7 ) 

6007 CONTINUE 

GO TO(610*611 ) *NPRT 

610 I CD= I CD+ 1 

CALL PLTSIX2 ♦ ICD»U2*V2*W2 *CP*RS2»MK) 

IF C ICD.GT.100) ICO*=0 

611 CONTINUE 

GO TO (2101 *9003 ), IFND 

9001 WRITE ( I W * 20022 ) 

20022 FORMAT ( 1 OX ♦ 42HSH0CK POINT CALCULATION NOT CONVERGING) 

GO TO 9999 

9002 WRITE (IW, 20023) K , J , DX 9 RP , VP2 ♦ A * EM ( J * K ) , BE TA , TANALF . ELAMB » TEST 

20023 FORMAT (///// 10X , 38HSECOND ORDER SOLUTION UNSTABLE K>l5*5X, 
12HJ= I5////9F13.8 ) 

GO TO 9999 

C PUNCH OUT RESULTS 

0003 GO TO ( 900A ,9999 ) , IPUNCH 

9004 WRITE (7,30001) < T < K ) , RS2 ( K ) , RSX2 ( K) , K = 1 ♦ MK > 

30001 FORMAT (6F12.7) 

DO 9005 K £ 1 « MK 

9005 WRITE (7,30002) ( U2 ( J ♦ K ) , V2 ( J *K ) , W2 ( J , X ) , P ( J ) , J = 1 , M J ) 

30002 FORMAT (4F15.7) 

WRITE(6»650)X2 

650 FORMAT ( *0 PUNCHED OUTPUT IS FOR X2 =*,E15.6) 

CALL PLTS(X2 *101 »U2»V2 *W2 »CP»RS2»MK) 

CALL GRPH(4,R2,R2,22J 
9999 CALL EXIT 
END 
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SUBROUTINE FIRST 

DIMENSION CONFIG! 17) *START{ 16) eT8(20) *KSP(20) ♦DTH(20) ,TD(20). 

1PHI (21*20) *TPHI (21*20) »FT(21.20)*U1(21»20) *V1 (21*20) » W1 ( 21 * 20 ) 
COMMON /FIR/XB ( 20 ) .RB(20*20> ,RBP( 20,20) »RBPP(20»20) »RBPPP ( 20* 20 ) * 
1UU( 20 ) » VU ( 20 ) » WU ( 20 ) ♦ Rl ( 20 ) ,RX 1 { 20 ) » RT 1 { 20 ) »DRl ( 20 ) *DRT1 (20) , 

2RS1 ( 20) »RSX1 ( 20) *RST1 (20) *F1(10080) *EL1(21*20»2)*T(20)»P(21)* 

3TOLFP.CVNST ,CTST .CNX.CNNX.XEND.DT.DTD.X1 » VI NF2 .G2.G3 .G4.G5 »CCP1 » 
4CCP2.MI .MJ.MK.MAXM.MAXN.MODE.IPT.KPT.MIBMl.MJMl 
COMMON /F2/S(50)/F12/FPPP(50)/F13/DS/F123/F(50) ,FP(50),FPP(50) 


1/C1/G1 ♦ AO/HESH/NOPT *ipunch 

EQUIVALENCE ( FI ( 1 ) »U1 ( 1 ) ) . ( FI ( 1681 ) « VI ( 1 ) ) » ( F 1 ( 3361 ) » W1 ( 1 ) ) * 
1(PHI(1)*TPHI(1)) 

IR = 5 
I W = 6 

READ ( I R * 10001 ) MIB.MKB»C0NFIG*?XB(TB)»IB*1»MIB) RED 

WRITE< IW. 10001 ) MIR. MK8 .CONFIG. <XB ( 18) ♦ 18=1 »MIB ) 

10001 FORMAT ( 2 I 5 » 17A4/ { 8F10 • 5 ) ) 

READ ( I R * 10002 ) ( TR( K8 ) *KB«1 »MKB ) RED 

WRITE( I W. 10002) ( T R ( K B ) ♦ KB'd.MKB) 

10002 FORMAT (8F10.5) 

DO 11 IB = 1 »MIB 

11 READ ( IR ♦ 10003 ) ( RB ( IB . KB ) .KB»1 .MKB ) RED 

10003 FORMAT (8F10.0) 

READ ( I R * 10004 ) ( RRP { 1 » KB ) . KB»1 .MKB ) RED 

10004 FORMAT (8F10o0) 

READ ( I R * 10005 ) ( RPP ( M I B * KB ) ♦ KB = 1 *MKB ) RED 

10005 FORMAT (8F10.0) 

READ (IR. 10006) MJD »MKD * I VAX I $ *ST ART » EMI NF , P I NF ,Rh0I NF , WHITE 

1 GAMMA* ALPHA. X1,(TD(KD> .RSl(KD) »RSX1 ( KD ) »KD« 1 *MKD ) WHITE 

WRITE( I W» 10006 ) MJD.MKD . I VAX I S .START , EMINF ,P INF »RHO I NF , WHITE 

1GAMMA. ALPHA. XI , (TD(KD) »RS1(KD) .RSXKKD) *KD* 1 »MKD ) WHITE 

10006 FORMAT ( 3 I 5 » 16A4/6F1 2 . 6 / ( 6F 1 2 0 6 ) ) 

DO 21 KD*1.MKD 

READ (IR. 10007) ( Ul ( JD ♦ KD ) » VI ( JD* KD ) » W1 < JD. KD ) . WHITE 

1 PH I ( JD * KD ) ♦ JD° 1 . M JD ) WHITE 

21 WRITE(IW. 10007) ( Ul ( JD . KD ) ♦ VI ( JD ♦ KD ) » W1 ( JD. KD ) * WHITE 

1PHI (JD.KD) »JD a l*MJD) 

10007 FORMAT (4F15.7) 

CALL NOMD I M ( Ul .VI ♦ W1 » GAMMA , EM I NF , P INF ♦ RHOINF »MJD »MKD ) 

READ (IR. 10008) MI ,MJ *MK ,MAXM .MAXN .MODE ♦ IPT » KPT * I PUNCH, NOPT , BLUE 

1T0LEP.CVNST .CTST .CNX.CNNX.XEND BLUE 

WR I TE ( I W» 10008 ) MI ,MJ ,MK ,MAXM .MAXN .MODE . IPT * KPT ♦ I PUNCH, NOPT , BLUE 

1T0LEP .CVNST .CTST .CNX.CNNX.XEND BLUE 

10008 FORMAT ( 10I5/6F10.5) 

WRITE (IW. 20001) CONFIG. START 


20001 FORMAT { 1H1 » 39X , 50HTHRFE DIMENSIONAL NEAR CHARACTERISTICS PROG 
1RAM///32X.13HCONFIGURATION, UX .17A4//32X.24HSTARTING VELOCITY DA 
2TA.16A4///) 

CONSTANT FACTORS (900) 

MIBM1*MIB-1 
MKBM1 *MK8-1 
MJDM1 *MJD-1 
MKDM1 “MKD-l 
MIM1»MI-1 
MJM1-MJ-1 
MKM1 *MK-1 
DP=1./MJM1 
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n r> 


P( 1 )=0.0 

DO 901 J=2«MJ 

901 P( J)=P< J-H J OP 
DTOR=3. 14159?/ 1=0 
DTD-TD(MKD) /MKMi 
DT = !'.TD*DTOR 

T { 1 ; 0 • 0 
DO 902 K = 2 »M< 

902 T(K)=T(K-1)+0TD 
EMTNF2=EMINF**2 

A I NF=6AMMA#P I NF / RH '• i 0 F 
AAINF = SQRT ( A INF ) 

VlN r 2=EMINF2*AINF 
VINF = S0RT(VINF2 ) 

61* • 5* ( GAMMA-1 <, J 

62= (GAMMA-1 . ) / ( GAMMA+l . 5 

G3 = 2.«-AINF/(GAMMA + 1 , ) 

G4= GAMMA/ ( GAMMA- 1 , ) 

G5 = 2 «*GAMMA/ ( A I NF* f 6AMMA+1 # ) ) 

VL I M2 = V I NF2+A I NF / G 1 
VLIM=SQRT(VLIM2 ) 

A0=G1*VLIM2 

CCP1=AINF**G4 

CCP2 c 2 . / ( GAMMA* EM INF2*CCP1 J 
IPUNCH-2-IPUNCH 

COMPLETE BODY RADIUS MATRIX (1000) 

IF (MKB-MK) 1001*1051*1001 

1001 DO 1002 KB=1»MK8 

1002 S ( KB ) =TB ( KB ) 

DO 1005 K=l,MKMl 
DO 1003 KB=1 » M K B M 1 
KSP ( K ) =KB 

IF (S(KB+1)-T(KIJ 1003.1003,1004 

1003 continue 

1004 KS =KSP ( K ) 

1005 DTH(K)«T(K)-S;kS) 

KSP (MK) =MKB 
DTH(MK) =0.0 

DO 1008 IB=1»MJB 
DO 1006 K8=l*MK3 

1006 F ( KB ) =RB ( IB. KB) 

CALL FIT2 (MKG.l ,1,0. -0. ) 

DO 1007 K = 1 ♦ MK 
KS=KSP(K) 

1007 RB( IB»K)=F(KS)+DTH(K)*(FP(KS)+.5*DTH(K)*(FPP(KS)+DTH(K)*FPPP(K5) / 
13. ) ) 

1008 CONTINUE 

DO 1009 KB= 1 , MKB 

1009 F ( KB ) =RBP ( 1 »K B ) 

CALL FIT2 (MKB.l ,1 .0. ,0. ) 

DO 1010 K=1,MK 
KS=KSP ( K ) 

1010 R0P(1»K)=F(KS) +DTH (K5#(FP(KS)+,5*DTH(K)*(FPP(KS) +DTH (K)*FPPP(KS)/ 
13. ) ) 

DO 1011 KB=1*MKB 

1011 F(KB)=RBP(MIB,K8) 

CALL FIT2 ( MKB *1*1 ,0. *0. ) 



DO 1012 K=1,MK " ’ 

KS=KSP(K) 

1°12 RBP (MIB»K ) =F ( KS)+DTH( K )*{ FP ( KS ) +, 5«DTH ( K ) * ( FPP ( KS ) +DTH ( K ) *FPPP (KS) 
1/3.)) 

1051 DO 1052 I B = 1 *M I B 

1052 S( I 8 ) =XB ( IB) 

DO 1055 K = 1 * MK 
DO 1053 I B= 1 ♦MIB 

1053 F ( I 8 ) =RB < I B ♦ K ) 

CALL f I T 2 (MI8,1*1,R8P(1*K) »R8P f Mf B »K ) ) 

DO 1054 I B= 1 *M I B 
RBP ( I B » K ) =FP ( IB) 

RBPP ( 1 8 » K ) = FPP { 18) 

1054 RBPPPf IB.K)=FPPP{ IB) , ■ 

10 55 CONTINUE ’-f 

WRITE ( I W, 20002) v.- 

20002 FORMAT ( ////40X » 30HCOMPLETE BODY RADIUS MATRIX/) 

DO 1056 IB=1»M!B 

1056 WRITE (IW, 20003) I B * XB ( I B ) * ( K *RB M B »K ) * RBP { I B ♦ K ) , RBPP ( I B » K ) , 

1 RBPPP ( 1 8 * K ) »K=1*MK) ri- ■■ 

20003 FORMAT (///10X»3HIR=I5*10X, 3HXB* F 10 . 5/ // 14X »lHK,15X,2HRB,l9x, 

1 3HRBP * 19X * 4HRBPP » 1 7X » 5HRBPPP /■/■( lOfi « I 5*4F22.6 ) ) 

C f . 

C UPSTREAM FLOW CONDITIONS ( 1200) 

VFACT1=VINF*C0S( ALPHA*DTOR) 

VFACT2=VINF*SIN( ALPHA*DTOR ) , 

DO 1201 K = 1 » MK 
UU ( K ) =VFACT 1 

VU(K)=-VFACT2*C0S(T(K)*DT0R) . 

1201 WU(K)*VFACT2*SINIT(K)*DT0R) 

WRITE (IW, 20004) EM I NF , P I NF ,RHOI NF ,G AMMA , AA I NF , VL I M , ALPHA , 

1(K»T(K) ,UU(K) *VU(K) »WU(K) *K*1*MK.) 

20004 FORMAT ( 1H1 , 40X , 26HUPSTREAM FLOW. COND I T I ONS/ // / 1 2X , 2 8HFREE STRE 
1AM MACH NUMBER *F 1 0 . 5 // 16X ♦ 24HFREE STREAM PRESSURE =F1C.5//17X 
2,23HFREE STREAM DENSITY *F 10. 8// 19X » 2 1HFREE STREAM GAMMa = 
3F10.5//8X»32HFREE STREAM SPEED -OF SOUND =F 1 0 . 4//20X » 20HL I M I T I N 
4G VELOCITY =F 10 . 4//2 IX , 1 9HANGLE OF ATTACK =F10.5////// 

519X , 1HK» 12X , 5HTHETA *20X ,2HUU*20X,2HVU, 20X »2HWU// ( 15X*I5»4F22.6) ) 

C 

C FIRST INITIAL DATA SURFACE (2000) 

DO 2002 18=1 ,MIBM1 

i s= i b . 

IF ( XB C IB + D-Xl) 2002*2002*2003 - 

2002 CONTINUE . 

2003 DEX = X1-XB( IS ) % 

DO 2004 K = 1 »MK ' 

R1 (K)=RB( IS,K)+DEX*( RBP( I5.K RBPP! I S » K ) +DEX *RBPPP ( IS,K) / 
13. ) ) 

RX1 ( K i =RBP( IS,K)+DEX# (RBPP( IS ,K >+.5*DEX*RBPPP ( IS*K ) ) 

2004 F(K)«R1 (K) 

DS = DT 

CALL F I T 3 ( MK * 1 ) 

DO 2005 K = 1 ,MK 

2005 RT1 ( K ) =FP ( K ) 

IF ( 1 - 1 VAX I S ) 2011*2020*2011 
2011 DO 2013 KD = 1 »MKD 
DO 2012 JD= 1 »MJD 
SPHI=SIN(PHI ( JD»KD)*DTOR) 
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2047 KEND=2 

2048 CALL FIT2 (MX 
DO 2049 K = l,f 
KS=KSP ( < ) 

2049 CALL VAL4(F1< 
1DTH(K) ) ) 

2050 CONTINUE 

2051 CONTINUE 

DO 2052 KD= 1 • 

2052 F(KD)=R51(KD: 
CALL FIT2 ( MK 
DO 2053 K=1*X 
KS=KSP ( K ) 

2053 RSl ( K ) =F ( KS H 
1 ) ) 

DO 2054 KD=1( 

2054 F ( KD ) =R$X 1 ( KC 
CALL FIT2 (MX 
DO 2055 K = 1,X 
KS=KSP ( K ) 

2055 RSX 1 ( K ) =F ( KS ) 
1 ) ) 

2061 DO 2062 K = 1*X 
DR 1 ( K ) = ( RSI ( X 

2062 F(K)=DR1(K) 
CALL FIT3 (MX 
DO 2064 K = 1 • N 
DRT 1 ( K ) *FP ( K ) 
RST 1 ( K ) =RT1 ( K 
DO 2063 JJ=1, 
CALL CHAR(U1( 

2063 CONTINUE 

2064 CONTINUE 
RETURN 
END 


CPH I =COS ( PH I ( JD* KD ) *DT OR ) 

TPHI ( JD*KD) =SPHI /COHI 

TEMP0R=U1 ( JD*KD)*SPHI+V1 ( JD»KD)*CPH 

U1 ( JD*KD)=U1 ( JD»<D)*CPHI-V1 ( JD*KD)* 

2012 VI ( JD»KD)=TEMPOR 

WR I T E ( I W* 10007 ) ( U 1 ( JD.KD) »V1 ( JD,KD 
1 PH I ( JD ♦ KD ) , JD=1 ,MJD) 

2013 CONTINUE 

2020 IF (MJD-MJ) 2021 , 2040 ,2021 

2021 DO 2032 KD=1 *VKD 
DO 2022 JD= 1 » M JD 

2022 S( JD) =TPHI ( J0*KD) 

DO 2025 NN = 1 *3 

DO 2023 JD= 1 » M JD 

2023 F(JD)=VAL(F1*1»JD*KD*1 *NN) 

CALL FIT2 (MJD,3,3,0. ,0. ) 

DO 2024 JD=1 *MJD 

CALL VAL4(F1 , 1 , JD *KD*2 »NN *FP ( JD) ) 
CALL VAL4(F1 ,1 , JD»KD»3 »NN*FPP{ JD) ) 

2024 CALL VAL4 ( F 1 * 1 , JD * KD * 4 *NN ♦ FPPP ( JD ) ) 

2025 CONTINUE 

DTPHIM SfMJD >-S(l ) ) /MJM1 
TP = S ( 1 ) -DTPHI 
DO 2029 J J= 1 »M J 

tp=tp+dtpht 

DO 2026 JD= 1 » M JDV 1 
JS = JD 

IF ( S ( JD+1 ) -TP ) 2026*2026*2027 

2026 CONTINUE 
JS=M JD 

2027 DTP=TP-S(JS) 

DO 2028 NN= 1,3 

2028 FT { J J *NN ) =HA RD I F ( VAL ( F 1 * 1 ♦ JS t KD * 1 t NN 
1 VAL ( F 1 * 1 »JS»KD»3»NN) * VAL ( F 1 * ] * JS ♦ KD » 

2029 CONTINUE 

DO 2031 J J= 1 * M J 
DO 2030 NN= 1*3 

2030 CALL VAL4(F1*1»JJ»KD»1*NN»FT(JJ»NN) ) 

2031 CONTINUE 

2032 CONTINUE 

2040 IF (MKD-MK) 2041*2061*2041 

2041 DO 2042 KD= 1 *MKD 

2042 S ( KD ) =TD ( KD ) 

DO 2045 K = 1 »MKMl 
DO 2043 KD= 1 ,MKDM1 
KSP (K )=KD 

IF ( S I KD+1 ) -T ( K ) ) 20^3*204^*2044 

2043 CONTINUE 

2044 KS=KSP(K) 

2045 DTH(K)=T ( K ) -S ( KS ) 

KSP ( MK ) =MKD 
DTH(MK)=0,0 

DO 2051 JJ=1 , M J 
DO 2050 NN= 1*3 
DO 2046 KD= 1 *MKD 

2046 F(KD)=VAL(F1 »1»JJ»KD*1*NN) 

KEND= 1 

IF (3-NN) 2047,2047*2048 
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FUNCTION HARDIF( A»B*C?D*DR) 

HARD! F* A+DR* ( B+O. 5*DR* < C+DR*D/3*0 ) ) 

RETURN 

END 
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FUNCTION VAL(VAR»IS*I » J ♦ K ♦ L ) 

DIMENSION VAR( 1 ) *NN( 3 *2 ) 

DATA NN/21»420*168O*21 *420*1260/ 

LOC=I + ( J-1)*NN( 1* I S > + C K — 1 } *NN ( 2 » I S ) + C L- 1 ) *NN ( 3 * I S ) 
IF(LOC.GT.10080-( IS-1 > * 3780 ) WR I TE ( 6 » 1001 ) ^ 

VAL =VAR ( LOC ) 

RETURN 

END 
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